\documentclass{article}
\usepackage{fullpage}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{courier}
\usepackage{appendix}


% If your system does not have the AMS fonts version 2.0 installed, then
% remove the useAMS option.
%
% useAMS allows you to obtain upright Greek characters.
% e.g. \umu, \upi etc.  See the section on "Upright Greek characters" in
% this guide for further information.
%
% If you are using AMS 2.0 fonts, bold math letters/symbols are available
% at a larger range of sizes for NFSS release 1 and 2 (using \boldmath or
% preferably \bmath).
%
% The usenatbib command allows the use of Patrick Daly's natbib.sty for
% cross-referencing.
%
% If you wish to typeset the paper in Times font (if you do not have the
% PostScript Type 1 Computer Modern fonts you will need to do this to get
% smoother fonts in a PDF file) then uncomment the next line
% \usepackage{Times}

%%%%% AUTHORS - PLACE YOUR OWN MACROS HERE %%%%%


%% cosmology
%% ions / absorbers
\newcommand{\HI}{\hbox{{\sc H}{\sc i}} }
\newcommand{\NHI}{{N_{\rm HI}}}
\newcommand{\fNHI}{f(N_{\rm HI},X)}

\newcommand{\DXobs}{\Delta X^{\rm obs}}
\newcommand{\DXsim}{\Delta X^{\rm sim}}

\newcommand{\dXobs}{dX^{\rm obs}}
\newcommand{\dXsim}{dX^{\rm sim}}

\newcommand{\nel}{n_{\rm e} }
\newcommand{\nH}{n_{\rm _{H}} }
\newcommand{\nHzero}{n_{\rm _{H,0}} }
\newcommand{\nHI}{n_{\rm _{HI}} }
\newcommand{\nHII}{n_{\rm _{HII}} }
\newcommand{\xHI}{x_{\rm _{HI}} }
\newcommand{\xHIzero}{x_{\rm _{HI,0}} }
\newcommand{\xHII}{x_{\rm _{HII}} }
\newcommand{\xHIIzero}{x_{\rm _{HII,0}} }
\newcommand{\xo}{x_{\rm _{0}} }

\newcommand{\sigth}{\sigma_{\rm 0}}
\newcommand{\Inu}{I_{\nu} }      % Specific Intensity
\newcommand{\nuth}{\nu_{\rm th}}
\newcommand{\numax}{\nu_{\rm max}}
\newcommand{\Jnu}{J_{\nu}}
\newcommand{\signu}{\sigma_{\nu}}
\newcommand{\taunu}{\tau_{\nu}}

\newcommand{\Ob}{\Omega_{\rm b} }
\newcommand{\Om}{\Omega_{\rm m} }
\newcommand{\OHI}{\Omega_{\rm HI} }
\newcommand{\OHH}{\Omega_{\rm H2} }
\newcommand{\Ol}{\Omega_{\Lambda} }
\newcommand{\Mpch}{h^{-1} \rm{\,Mpc} }

% codes / simulations
\newcommand{\owls}{{\small OWLS} }
\newcommand{\urchin}{{\small URCHIN} }

% journals

\newcommand{\mnras}{MNRAS}
\newcommand{\apjl}{ApJL}
\newcommand{\apj}{ApJ}
\newcommand{\apjs}{ApJS}
\newcommand{\aap}{A\&A}
\newcommand{\araa}{ARA\&A}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title{Hydrogen Ionization States \\ Analytic Solutions}
\author{Gabriel Altay}
\begin{document}
\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

These notes are meant to document a series of analytic solutions for the evolution of the ionization state of hydrogen.  They consider a simple model atom that can only be in the ground state or ionized. 








%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Basic Equations}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The evolution of the hydrogen neutral fraction, $\xHI = \nHI/\nH$, or equivalently ionized fraction, $\xHII = \nHII/\nH$, is a balance between photo/collisional ionizations and recombinations.

\begin{eqnarray}
\label{dx1dt}
\frac{d\xHI}{dt} &=& -(\Gamma + \gamma \nel) \xHI + 
\alpha  \nel (1 - \xHI) \\
\label{dx2dt}
\frac{d\xHII}{dt} &=& (\Gamma + \gamma \nel)(1 - \xHII) - 
\alpha \nel \xHII
\end{eqnarray}
where $\Gamma$ is the photoionization rate, $\gamma$ is the collisional ionization rate, $\alpha$ is the recombination rate, and $\nel$ is the free electron number density.  If we decompose the free electron number density into a contribution from ionized hydrogen and a contribution from all heavier elements, we can write $\nel = (\xHII + y) \nH$.


\section{Equilibrium - Pure Hydrogen - Cold Gas}

In this approximation $y=0$, $\nel = \xHII \nH = (1-\xHI) \nH$, $\gamma = 0$ and Eqs.~\ref{dx1dt} and \ref{dx2dt} reduce to,
\begin{eqnarray}
0 &=& -\Gamma \xHI + 
\alpha \nH (1-\xHI)^2 \\
0 &=& \Gamma (1 - \xHII) - 
\alpha \nH \xHII^2
\end{eqnarray}
\\
\subsection{Neutral Fraction}
Grouping terms in powers of $\xHI$ and defining $U \equiv \Gamma / (\nH \alpha)$ we have,
\begin{eqnarray}
0 &=& -\Gamma \xHI + \alpha \nH (1 - \xHI)^2 
\\
  &=& -U \xHI + (1 - 2\xHI + \xHI^2)
\\
  &=& \xHI^2 - (U + 2) \xHI + 1
\end{eqnarray}
\\
The physical solution takes the negative sign in the quadratic equation
\begin{eqnarray}
2 \xHI &=& (U+2) - \left[ (U+2)^2 - 4 \right]^{1/2} \\
2 \xHI &=& (U+2) - \left[ U^2 + 4U \right]^{1/2} 
\end{eqnarray}
\\
In a highly neutral gas $\Gamma \ll \alpha \nH$ and $U \ll 1$. Using the binomial approximation we have, 
\begin{eqnarray}
2 \xHI &=& (U+2) - 2U^{1/2} \left[1 + \frac{1}{4} U \right]^{1/2} \\
\xHI &=& \left( \frac{1}{2}U+1 \right) - U^{1/2} 
\left[1 + \frac{1}{8} U - \frac{1}{128} U^{2} + ...\right] \\
\xHI &=& 1 - U^{1/2} + 
\frac{1}{2} U - \frac{1}{8} U^{3/2} + \frac{1}{128} U^{5/2} + \mathcal{O}(U^{7/2})
\end{eqnarray}
\\
In a highly ionized gas $\Gamma \gg \alpha \nH$ and $U^{-1} \ll 1$. Using the binomial approximation  we have, 
\begin{eqnarray}
2 \xHI &=& (U+2) - U \left[1 + 4U^{-1} \right]^{1/2} \\
2 \xHI &=& (U+2) - U \left[1 + 2 U^{-1} - 2 U^{-2} + ...\right] \\
2 \xHI &\approx&  U+2 - U - 2 + 2 U^{-1} \\
\xHI &\approx& U^{-1} = \frac{\alpha \nH}{\Gamma} + \mathcal{O}(U^{-2})
\end{eqnarray}



\section{Exact Solutions}

Eqs. 1 and 2 can both be put into the form of a Riccati equation, 
\begin{eqnarray}
\frac{dx}{dt} = R x^2 + Q x + P.
\end{eqnarray}
Below we will determine the coefficients $R, Q, $ and $P$ for the neutral and ionized fractions. 


\subsection{Neutral Fraction Coefficients}

Grouping the terms in powers of $\xHI$, the first term in Eq. 1 becomes,
\begin{eqnarray} 
-(\Gamma + \gamma \nel) \xHI = \\
-\Gamma \xHI - \gamma [(1-\xHI) + y] \nH \xHI = \\
-\Gamma \xHI -\gamma \nH \xHI + \gamma \nH \xHI^2 - \gamma \nH y \xHI = \\
\gamma \nH \xHI^2  
- (\Gamma + \gamma \nH + \gamma \nH y ) \xHI
\end{eqnarray}
and the second becomes,

\begin{eqnarray} 
\alpha \nel (1 - \xHI) = \\
\alpha [(1-\xHI) + y] \nH (1-\xHI) = \\
\alpha \nH (1-\xHI)^2 + \alpha \nH y (1-\xHI) = \\
\alpha \nH (1 - 2\xHI + \xHI^2) + \alpha \nH y - \alpha \nH y \xHI =\\
\alpha \nH \xHI^2 - (2 \alpha \nH + \alpha \nH y) \xHI + 
\alpha \nH (1 +y)
\end{eqnarray}
combining the two we get, 

\begin{eqnarray}
(\gamma + \alpha) \nH \xHI^2 - 
\left[ \Gamma + \gamma \nH + \gamma \nH y  
+ 2 \alpha \nH + \alpha \nH y \right]  \xHI + 
\alpha \nH (1+y) = 0 \\
(\gamma + \alpha) \nH \xHI^2 -
\left[ \Gamma + (\gamma + 2\alpha) \nH + 
( \gamma + \alpha ) \nH y 
\right] \xHI + 
\alpha \nH (1+y) = 0 
\end{eqnarray}
Matching coefficients we get, 

\begin{eqnarray}
R &=& (\gamma + \alpha) \nH  \\
Q &=& -\left[ \Gamma + (\gamma + 2\alpha) \nH + 
(\gamma + \alpha) \nH y \right]  \\
  &=& -\left[ \Gamma + \alpha \nH + R (1+y) \right]  \\
P &=& \alpha \nH (1+y)  
\end{eqnarray}
Note that in the case of the neutral fraction, all the coefficients have a definite sign no matter what value the parameters take.  


\subsection{Ionized Fraction Coefficients}


Grouping the terms in powers of $\xHII$, the first term in Eq. 2 becomes,

\begin{eqnarray} 
(\Gamma + \gamma \nel) (1-\xHII) = \\
\left[ \Gamma + \gamma (\xHII + y) \nH \right]  (1-\xHII)  = \\
\Gamma + \gamma \nH \xHI + \gamma \nH y  -
\Gamma \xHI - \gamma \nH \xHI^2 - \gamma \nH y \xHI = \\
- \gamma \nH \xHI^2 + 
( \gamma \nH - \gamma \nH y ) \xHI
+ \Gamma + \gamma \nH y  
\end{eqnarray}
and the second becomes

\begin{eqnarray} 
- \alpha \nel \xHII = \\
- \alpha (\xHII + y) \nH \xHII = \\
-\alpha \nH \xHII^2 - \alpha \nH y \xHII
\end{eqnarray}
combining the two we get, 

\begin{eqnarray}
-( \gamma \nH + \alpha \nH ) \xHII^2 + 
( \gamma \nH - \Gamma - \gamma \nH y - \alpha \nH y) \xHII
+ \Gamma + \gamma \nH y = \\
-( \gamma + \alpha ) \nH \xHII^2 + 
\left[ \gamma \nH - \Gamma - (\gamma + \alpha) \nH y \right] \xHII
+ \Gamma + \gamma \nH y 
\end{eqnarray} 
Matching coefficients we get, 

\begin{eqnarray}
  R &=& -(\gamma + \alpha) \, \nH \\
  Q &=& \gamma \nH - \Gamma - 
  (\gamma + \alpha) \, \nH \, y  \nonumber \\
    &=& \gamma \nH - \Gamma + R \, y  \\
  P &=& \Gamma + \gamma \, \nH \, y   
\end{eqnarray}
In this case, the sign of $Q$ depends on the choice of parameters. 



\subsection{Solutions}

Now we have defined the coefficients $R$, $Q$, and $P$, we will explore solutions for $\xHI$ and $\xHII$.  Note that the following applies to both solutions (neutral and ionized fraction).  The only thing that changes is the values $R$, $Q$, and $P$.   

\subsubsection{Approach 1}
We first note that any quadratic expression can be factored in terms of its roots $x_{+}$ and $x_{-}$
\begin{eqnarray}
\frac{dx}{dt} = R x^2 + Q x + P = R ( x - x_{+} ) ( x - x_{-} ) 
\end{eqnarray}
where, 
\begin{eqnarray}
x_{\pm} = \frac{-Q \pm (Q^2 - 4 R P)^{1/2}}{2R}, \quad {\rm and} \quad
\delta x_\pm \equiv x_{+} - x_{-} = 
\frac{ (Q^2 - 4 R P)^{1/2}}{R}
\end{eqnarray}
The physical equilibrium solution is $x_{\rm eq} = x_{-}$.  To find the time dependent solution we integrate.

\begin{eqnarray}
\int_{x_0}^x \frac{dx'}{R ( x' - x_{+} ) ( x' - x_{-} ) } &=& \int_0^t dt 
\end{eqnarray}

\subsubsection{Method 1 - Partial Fractions}

\begin{eqnarray}
\frac{1}{(x-x_{+})(x-x_{\rm eq})} &=& 
\frac{A}{x-x_{+}} + \frac{B}{x-x_{\rm eq}} 
\\
1 &=&
A (x-x_{\rm eq}) + B (x-x_{+})
\end{eqnarray}
These relations must hold for any value of $x$.

\begin{eqnarray}
{\rm for} \quad x = x_{+}, &&\quad A = \frac{1}{x_{+}-x_{\rm eq}} 
= \frac{1}{\delta x_\pm} 
\\
{\rm for} \quad x = x_{\rm eq}, &&\quad B = \frac{1}{x_{\rm eq}-x_{+}} = -A
\end{eqnarray}
therefore, 

\begin{eqnarray}
  \frac{1}{(x-x_{+})(x-x_{\rm eq})} = 
  A \left( \frac{1}{x-x_{+}} - \frac{1}{x-x_{\rm eq}} \right)
\end{eqnarray}

\begin{eqnarray}
  \int_{x_0}^x \frac{dx'}{R ( x' - x_{+} ) ( x' - x_{-} ) } &=& \int_0^t dt 
  \\
  \frac{A}{R} \left( 
    \int_{x_0}^x \frac{dx'}{x' - x_{+}} - 
    \int_{x_0}^x \frac{dx'}{x' - x_{\rm eq}} 
  \right)  &=& t
  \\
  \left. \ln(x-x_{+}) \right|_{x_0}^x -  
  \left. \ln(x-x_{\rm eq}) \right|_{x_0}^x -  
   &=& \frac{R t}{A}
   \\
   \ln \left( \frac{ x-x_{+} }{ x_0 - x_{+} } \right) - 
   \ln \left( \frac{ x-x_{\rm eq} }{ x_0 - x_{\rm eq} } \right) 
   &=& R \delta x_{\pm} t
   \\
   \ln \left[ \left(  \frac{ x-x_{+} }{ x_0 - x_{+} } \right) 
   \left( \frac{ x_0 - x_{\rm eq} }{ x-x_{\rm eq} } \right) \right] 
   &=& R \delta x_{\pm} t
   \\
   \left(  \frac{ x-x_{+} }{ x_0 - x_{+} } \right) 
   \left( \frac{ x_0 - x_{\rm eq} }{ x-x_{\rm eq} } \right) 
   &=& \exp( R \delta x_{\pm} t )
   \\
   \left( \frac{x - x_{\rm eq}}{x - x_{+}} \right)    
   \left( \frac{x_0 - x_{+}}{x_0 - x_{\rm eq}}  \right)
   &=& \exp \left( -R \delta x_{\pm} t \right)
   \\
   \left( x - x_{\rm eq} \right) \left(x_0 - x_{+} \right)      
   &=& \left( x - x_{+} \right) \left( x_0 - x_{\rm eq}\right) 
   \exp \left( -R \delta x_{\pm} t \right)
\end{eqnarray}

Defining some helpful notation,
\begin{eqnarray}
  t_{\rm eq} &=& (R \delta x_{\pm})^{-1} = (Q^2 - 4 R P)^{-1/2} \\
  \mathcal{F} &=& \exp \left( - t/t_{\rm eq} \right)
\end{eqnarray}

\begin{eqnarray}
  \left( x - x_{\rm eq} \right) \left(x_0 - x_{+} \right)      
  &=& \left( x - x_{+} \right) \left( x_0 - x_{\rm eq}\right) 
  \mathcal{F}
  \\
  x        \left(x_0 - x_{+} \right) - 
  x_{\rm eq} \left(x_0 - x_{+} \right)      
  &=& x    \left( x_0 - x_{\rm eq}\right) \mathcal{F} - 
  x_{+} \left( x_0 - x_{\rm eq}\right) \mathcal{F} 
  \\
  x \left[ \left(x_0 - x_{+} \right) - 
    \left( x_0 - x_{\rm eq}\right) \mathcal{F} \right]   
  &=& x_{\rm eq} \left(x_0 - x_{+} \right) -       
  x_{+} \left( x_0 - x_{\rm eq}\right) \mathcal{F} 
  \\
  x \left[ \left(x_0 - x_{+} \right) - 
    \left( x_0 - x_{\rm eq}\right) \mathcal{F} \right]   
  &=& x_{\rm eq} \left[ \left( x_0 - x_{+} \right) - 
    \left( x_0 - x_{\rm eq} \right) \mathcal{F} \right] + 
  \\
  && 
  x_{\rm eq} \left( x_0 - x_{\rm eq} \right) \mathcal{F} -       
  x_{+} \left( x_0 - x_{\rm eq}\right) \mathcal{F} 
  \\
  x &=& x_{\rm eq} +  
  \frac{ 
    x_{\rm eq} \left( x_0 - x_{\rm eq} \right) \mathcal{F} -       
    x_{+} \left( x_0 - x_{\rm eq}\right) \mathcal{F} }
  {
    \left(x_0 - x_{+} \right) - 
    \left( x_0 - x_{\rm eq}\right) \mathcal{F} }
  \\
  x &=& x_{\rm eq} +  
  \frac{ 
    \left( x_{\rm eq} - x_{+} \right) 
    \left( x_0 - x_{\rm eq} \right) \mathcal{F}  }
  {
   x_0 (1 - \mathcal{F} ) + \mathcal{F} x_{\rm eq} - x_{+} }
  \\
  x &=& x_{\rm eq} +  \left( x_0 - x_{\rm eq} \right)
  \frac{ 
    \left( x_{\rm eq} - x_{+} \right) 
     \mathcal{F}  }
  {
   x_0 (1 - \mathcal{F} ) + \mathcal{F} x_{\rm eq} - x_{+} }
\end{eqnarray}

Another approach is the following, 


\begin{eqnarray}
  \int_{x_0}^x \frac{dx'}{R ( x' - x_{+} ) ( x' - x_{-} ) } &=& \int_0^t dt 
  \\
  \left.  
    \frac{\log(x'- x_{+}) - \log(x'-x_{-}) }{\delta x_{\pm}}
  \right|_{x_0}^x 
  &=& R t 
  \\ 
  \left. \log \left( \frac{x'-x_{+}}{x'- x_{-}}  \right) \right|_{x_0}^x 
  &=& R \delta x_{\pm} t 
  \\ 
  \log \left( \frac{x   - x_{+}}{x   - x_{-}}  \right) -
  \log \left( \frac{x_0 - x_{+}}{x_0 - x_{-}}  \right)
  &=& R \delta x_{\pm} t 
  \\ 
  \log \left[ 
    \left( \frac{x   - x_{+}}{x   - x_{-}} \right)    
    \left( \frac{x_0 - x_{-}}{x_0 - x_{+}}  \right)
  \right]
  &=& R \delta x_{\pm} t 
  \\ 
  \left( \frac{x   - x_{+}}{x   - x_{-}} \right)    
  \left( \frac{x_0 - x_{-}}{x_0 - x_{+}}  \right)
  &=& \exp \left( R \delta x_{\pm} t \right)
\end{eqnarray}
Taking the reciprocal and rearranging we get: 
\begin{eqnarray}
\left( \frac{x   - x_{-}}{x   - x_{+}} \right)    
\left( \frac{x_0 - x_{+}}{x_0 - x_{-}}  \right)
 &=& \exp \left( -R \delta x_{\pm} t \right)
\\ 
\left( \frac{x   - x_{-}}{x   - x_{+}} \right)    
 &=& \left( \frac{x_0 - x_{-}}{x_0 - x_{+}}  \right)
\exp \left( -R \delta x_{\pm} t \right) 
\end{eqnarray} 
Replacing $x_{-}$ with $x_{\rm eq}$ and defining some helpful notation,
\begin{eqnarray}
\mathcal{D} &=& \frac{x_0 - x_{\rm eq}}{x_0 - x_{+}} \\
t_{\rm eq} &=& (R \delta x_{\pm})^{-1} = (Q^2 - 4 R P)^{-1/2} \\
\mathcal{F} &=& \exp \left( - t/t_{\rm eq} \right)
\end{eqnarray}
we can write

\begin{eqnarray}
\left( \frac{x - x_{\rm eq}}{x - x_{+}} \right) &=& \mathcal{D} \mathcal{F}
\\
x - x_{\rm eq} &=& \left( x - x_{+} \right) \mathcal{D} \mathcal{F}
\\ 
x - x \mathcal{D} \mathcal{F} &=& 
x_{\rm eq} - x_{+} \mathcal{D} \mathcal{F}
\\ 
x \left( 1 - \mathcal{D} \mathcal{F} \right) &=& 
x_{\rm eq} - x_{+} \mathcal{D} \mathcal{F}
\\ 
x(t) &=& 
\frac{ x_{\rm eq} - x_{+} \mathcal{D} \mathcal{F}(t) }
{1 - \mathcal{D} \mathcal{F}(t) }
\end{eqnarray} 
This form makes it clear that as $t \rightarrow \infty$ and $\mathcal{F}(t) \rightarrow 0$, the solution $x(t) \rightarrow x_{\rm eq}$.  If we expand $\mathcal{D}$ we can present a form that makes it clear that as $t \rightarrow 0$ and $\mathcal{F}(t) \rightarrow 1$, the solution $x(t) \rightarrow x_0$.
\begin{eqnarray}
x &=& 
\frac{ x_{\rm eq}  - x_{+} \frac{x_0 - x_{\rm eq}}{x_0 - x_{+}} \mathcal{F} }
{1 - \frac{x_0 - x_{\rm eq}}{x_0 - x_{+}} \mathcal{F} }
\\
x &=& 
\frac{ x_{\rm eq} (x_0 - x_{+}) - x_{+} (x_0 - x_{\rm eq}) \mathcal{F} }
{(x_0 - x_{+}) - (x_0 - x_{\rm eq})  \mathcal{F} }
\\
x &=& 
\frac{ 
x_{\rm eq} x_0 - x_{\rm eq} x_{+} - x_{+} x_0 \mathcal{F} + x_{+} x_{\rm eq} \mathcal{F} }
{x_0 - x_{+} - x_0 \mathcal{F} + x_{\rm eq} \mathcal{F} }
\\
x(t) &=& 
\frac{ 
x_0 \left[ x_{\rm eq}  - x_{+} \mathcal{F}(t) \right]
- x_{\rm eq} x_{+} [1-\mathcal{F}(t)  ] }
{ [ x_{\rm eq} \mathcal{F}(t) - x_{+}] + x_0 [1- \mathcal{F}(t) ]  }
\end{eqnarray}

\subsubsection{Approach 2}

The time dependent solution of this Ricatti equation is, 
\begin{eqnarray}
 x(t) = \frac{1}{2R} \left[    
   -Q - d \frac{\tanh(t/t_{\rm c}) - x_{\rm s}}
   {1 - x_{\rm s} \tanh(t/t_{\rm c})}
   \right] \nonumber \\
  t_{\rm c} = \frac{2}{d}, \quad 
  x_{\rm s} = \frac{Q + 2 R \xo}{d}, \quad
  \xo = x(0), \quad
  d = \sqrt{Q^2 - 4 R P}
\end{eqnarray}
and the physical equilibrium solution is given by the negative root of
the quadratic formula.  

\begin{equation}
x^{\rm eq} = \frac{-Q -\sqrt{Q^2 - 4 R P}}{2R}
\end{equation}
In a gas composed only of Hydrogen and Helium, $y$ is bound between
zero and (1-$X$)/(2$X$) where $X$ is the Hydrogen mass fraction of the
gas.  For this solution then the ionization state is a function of
four variables $\nH$, $T$, $\Gamma$, and $y$.  This solution is best
evaluated numerically using the following form, 
\begin{equation}
x^{\rm eq} = \frac{P}{q}, \quad {\rm with} \quad
q = -\frac{1}{2} \left[ Q + {\rm sgn}(Q) \sqrt{Q^2 - 4RP} \right]
\end{equation}

\subsubsection{Physical Intuition}

We can gain some physical intuition about the parameters by considering some simplified cases.  Examining the neutral fraction has the great benefit that the coefficients $R$, $Q$, and $P$ have fixed signs.  In a pure hydrogen gas ($y=0$) with with no incident radiation ($\Gamma=0$) we have (for the neutral fraction case), 

\begin{eqnarray}
P = \alpha \nH  > 0, \quad 
R = P + \gamma \nH > P, \quad 
-Q = R + P
\end{eqnarray}
which makes it clear that $P < R < -Q$.  It follows that 
\begin{eqnarray}
Q^2 &=& (R+P)^2 \\
Q^2 &=& R^2 + P^2 + 2RP \\
Q^2 - 4RP &=& R^2 + P^2 - 2RP \\
Q^2 - 4RP &=& (R - P)^2
\end{eqnarray}
As $(R-P)^2$ must be positive, this implies that $Q^2 > 4RP$ and that $x_{+} > x_{-}$.  If we consider the case with a finite photoionization rate $\Gamma$ then the only coefficient that changes is $Q \rightarrow Q + \Gamma$.  If $Q^2 > 4RP$ is true then so is $(Q+\Gamma)^2 > 4RP$. 



\appendix
\appendixpage
\addappheadtotoc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Basic Math}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Simple transformations between linear and log variables. 

\begin{equation}
\frac{d\ln x}{dx} = 
\frac{1}{x}, \quad
d\ln x = \frac{dx}{x}
\end{equation}

\begin{equation}
\frac{d\log x}{dx} = 
\frac{1}{\ln 10} \frac{d\ln x}{dx} =
\frac{1}{\ln 10} \frac{1}{x}, \quad
d\log x = \frac{1}{\ln 10} \frac{dx}{x}
\end{equation}
\\
\\
Binomial approximation for small $x$, 
\begin{equation}
(1+x)^n \approx 1 + n x + \frac{1}{2} n (n-1) x^2 + ...
\end{equation}


\begin{figure}
  \begin{center}
    \includegraphics[width=0.45\textwidth]{python/gamma_over_alpha}
  \end{center}
  \caption{Ratio of hydrogen collisional ionization rate ($\gamma$) and recombination rate ($\alpha$).  The case A recombination rate is shown in blue and case B is shown in red.  The ratio covers approximately eight orders of magnitude between $T = 10^4$ K and $T=10^5$ K. We have used the fits from \cite{1997MNRAS.292...27H}.  }
\end{figure}




\bibliographystyle{alpha} % or "unsrt", "alpha", "abbrv", etc.
\bibliography{biblio}	  % use data in file "astrobibl.bib"


\end{document}
